clear;
close all;
warning('off','all');
paths;

load('initial.mat');
load('data.mat');

T = height(data);

infl_r = NaN(T,1); alpha_r = NaN(T,1); phi_r = NaN(T,1);
xg_r = NaN(length(x0),T); 
ms = MultiStart;

for sss = 20:T

    display(data.YYYYQ(sss));

    problem = createOptimProblem('fminunc','x0',x0,'objective',@(x)obj_main(x, psi0, PP0, y(1:sss,:), y5(1:sss,:), pi(1:sss), 0),'options',optimoptions(@fminunc,'Algorithm','quasi-newton','Display','off','MaxFunctionEvaluations',50000,'MaxIterations',5000));
    [xg_r(:,sss), ~] = run(ms,problem,50);
    [~,est] = measurement(xg_r(1:4, sss), xg_r(5:end, sss), psi0, PP0, y(1:sss,:), y5(1:sss,:), pi(1:sss), 0);
    infl_r(sss) = est.shat(1,1,end);
    alpha_r(sss) = est.shat(2,1,end);
    phi_r(sss) = est.shat(3,1,end); 

end

phi_r = 1./(1+0.1.*exp(-phi_r));
realtime = alpha_r + phi_r .* infl_r;
soln.mse = mean((realtime(20:end-1) - pi(21:end)).^2);
soln.mae = mean(abs(realtime(20:end-1) - pi(21:end)));